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A general approach is presented for understanding the stress response function in anisotropic 
granular layers in two dimensions. The formalism accommodates both classical anisotropic elasticity 
theory and linear theories of anisotropic directed force chain networks. Perhaps surprisingly, two- 
peak response functions can occur even for classical, anisotropic elastic materials, such as triangular 
networks of springs with different stiffnesses. In such cases, the peak widths grow linearly with 
the height of the layer, contrary to the diffusive spreading found in 'stress-only' hyperbolic models. 
In principle, directed force chain networks can exhibit the two-peak, diffusively spreading response 
function of hyperbolic models, but all models in a particular class studied here are found to be in 
the elliptic regime. 

PACS numbers: 45.70.Cc Static sandpiles; granular compaction - 83.80.Fg Granular solids 



I. INTRODUCTION 

The stress response of an assembly of hard, cohesionlcss grains has been a subject of debate 
01 Q ■ The dividing line has been mostly between traditional approaches based on elasticity or 
elastoplasticity theory on one hand, and "stress-only" models on the other which make no reference 
to a local deformation field but posit (history-dependent) closure relations between components 
of the stress tensor. The former leads to elliptic partial differential equations for the stresses, for 
which boundary conditions must be imposed everywhere on the boundary. In contrast, the latter 
approach often leads to hyperbolic equations 0, Q . The wave- like behavior of their solutions has 
been at the origin of a proposed physical mechanism called stress propagation through the bulk 
a granular material. In an infinite slab geometry, it only requires the specification of boundary 
conditions on the "top" surface. A family of (linear) closure relations have been shown to account 
for the pressure dip underneath the apex of a sandpile and stresses in silos 0, ■ Alternative 
explanations based on elastoplasticity are found in pj. 

The phenomenological 'stress-only' closure relations follow from plausible symmetry arguments, 
and can be seen as the coarse-grained version of local probabilistic rules for stress transfer [|j . How- 
ever, these relations lack a detailed microscopic derivation that would allow one both to understand 
their range of validity and to compute the phenomenological parameters from the statistical prop- 
erties of the packing, except in the case of frictionle ss g rains. In fact, a system of frictionless 
polydisperse spheres is shown to be isostatic 0, H, El 1 1 0l | . i.e. the number of unknown forces is 
equal to the number of equations for mechanical equilibrium. If an isostatic system is sufficiently 
anisotropic, a linear closure relation between stresses can be derived [TlT |. Further attempts to 
obtain the missi ng e quation for stresses from a microscopic approach for different packings are 
presented in 0, [Q, but these are still somewhat inconclusive. In particular, in the case of a 
completely isotropic packing, none of the homogeneous linear closure relations is compatible with 
the rotational symmetry. The idea of 'grains' (in the metallurgical sense) and packing defects must 
be introduced to restore the large scale symmetry. 

In order to understand stress distribution on a more fundamental level, we have introduced the 
mesoscopic concept of the directed force chain network (dfcn) 0,0], which is motivated by the 
experimental evidence for filamentary force chains in a wide variety of systems. [l^. The "double 
Y" model has been developed to describe such networks based on simple rules for the splitting 
and merging of straight force chains. This model leads to a non-linear Boltzmann equation for the 
probability P(f, n, r) of finding a force chain at the spatial point r with intensity / in the direction 
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In a first paper [Tjj, chain merging (which produces the non- linear terms in the Boltzmann 
equation) was neglected. An isotropic splitting rule was assumed, corresponding to strongly dis- 
ordered isotropic granular packings. A pseudo-elastic theory for the stress tensor was derived in 
which the role of the displacement field is played by a vector field J(r) = (n/) that represents the 
coarse-grained or ensemble averaged force chain direction. A relation between di Jj and the stress 
tensor exists that is formally equivalent to an isotropic stress-strain relation. The resulting elliptic 
equations yield a response function with a unique (pseudo-elastic) peak, as observed experimen- 
tally in strongly disordered packings |l6lll7l|. Further study showed, however, that the non-linear 
terms in the Boltzmann equation contain essential physics and cannot be neglected. 01 I n fact, 
for an exactly solvable model with 6 discrete directions for force chains, it was found that the 
elliptic (pseudo-elastic) behavior of the response function is limited to small depths, and that at 
sufficiently large depths a crossover occurs to an hyperbolic response - two Gaussian peaks that 
propagate away and broaden diffusively. Whether this behavior is specific to the model with 6 dis- 
crete directions is a subject of current study, and the elliptic or hyperbolic nature of the linearized 
response around the full solution of the nonlinear Boltzmann equation is an open question. 

Following a different route, Goldcnberg and Goldhirsch [l^ have recently noted that a two-peak 
resp onse function can be found in classical anisotropic ball-and-spring models. Gay and da Silveira 
|l9j have furthermore given some arguments for the relevance of anisotropic elasticity for the large 
scale description of granular assemblies of compressible grains that can locally rotate. The two- 
peak nature of the response function is therefore not, in itself a signature of hyperbolicity, but may 
occur in elliptic systems that are sufficiently anisotropic. The unambiguous signature of hyperbolic 
response lies in the scaling of the peak widths with depth, which is linear in generic elliptic systems 
but diffusive (proportional to the square root of depth) in generic hyperbolic systems. In the linear 
pseudo-elasticity theories discussed below, the diffusive spreading in hyperbolic systems is not 
captured; the peaks appear as delta functions that do not spread at all. Deviations from elasticity 
on small scales and their possible relation with granular media were also discussed in |2(il | . 

The aim of this paper is to give a unified account of the shape of the response function for 
anisotropic systems described either by standard elasticity theory or the pseudo-elastic theory that 
emerges from an approximate linear treatment of directed force networks. Though there are open 
questions concerning the self-consistency of the latter, there do appear to be some contexts in 
which the equations of the pseudo-elasticity theory hold, and they may be especially relevant for 
systems of intermediate depth (large compared to the disorder length scale but not much larger 
than the persistence length of force chains). 

Very recently, the response functions of two-dimensional granular layers subjected to shear have 
been determined experimentally [27]]. Under shear, an anisotropic texture appears and force chains 
are preferably oriented along an angle of 45 degrees. Within the (pseudo)-elasticity framework 
presented below, this provides motivation for studying materials characterized by a selected global 
direction N. 

The paper is organized as follows. In section 2 a general mathematical framework for calculating 
stress response functions in anisotropic materials. The main results of the paper are then summa- 
rized in a 'phase-diagram' indicating where 'one-peak' or 'two- peak' response functions can appear 
in parameter space. In section 3, we compute the analytic form of the response function for the 
various phases and show a number of examples of the variety of shapes that are possible, includ- 
ing a brief comment on relation to experimental work. In section 4, we show how the formalism 
applies to the example of a triangular ball-and-spring network, indicating how spring stiffnesses 
must be chosen to access all possible regions of the general parameter space. In section 5, a linear 
anisotropic pseudo-elastic theory is derived from an anisotropic linear directed force chain network 
model and it is shown that this class of models always lies in the elliptic regime. A conclusion is 
given in section 6. Algebraic details of several calculations are presented in Appendices. 

II. ANISOTROPIC ELASTICITY AND SUMMARY OF OUR RESULTS 

A. General equations for 2D systems with arbitrary 
anisotropy 

In the following, we present a general framework that covers both classical linear anisotropic 
elasticity theory |2l| and a generally anisotropic "pseudo-elasticity" theory, that appears within 
a linearized treatment of directed force chain networks (see section 01 . The large scale equations 
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that can be derived in these two approaches are formally identical, although the "pseudo-strain" 
has a geometric meaning different from the usual strain tensor. For simplicity, we will restrict the 
discussion to two-dimensional systems. 

The most general linear relation between the stress tensor a and a symmetric tensor formed 
from the gradients of a vector field u is 

0~ij = KjklUkU (l) 

where ay denotes a component of the stress tensor, uy = ^(djUi + diUj), and summation over 
repeated indices is implied. In the classical linear theory of elasticity, the vector u, is the displace- 
ment field describing the physical deformation of a continuous medium. For usual elastic bodies, 
the antisymmetric combination djUi — diUj corresponds to a local rotation of the material, which is 
not allowed here. For granular materials, on the other hand, grains might locally rotate due to the 
presence of friction. This extension which suggests a continuum description in terms of Cosserat 
elasticity was recently discussed in |T^ . The absence of internal torques requires that the stress 
tensor is also symmetric. The coefficients Xijki are material constants and form the elastic modulus 
tensor. The indices i,j, k, I are equal to x, z where for later purposes x is to be considered as the 
horizontal coordinate and z a vertical coordinate pointing downward. 

Symmetry of both the stress and the strain tensor imply a permutation symmetry within the 
first and second pair of indices for Xijki, i.e. 

Kjki = Xjiki = Xijik = Xjiik- (2) 

Materials whose behavior is modeled only in terms of equation without any reference to a free 
energy functional are characterized by an elastic modulus tensor that need not have any symmetries 
other than equation They arc called 'hypoclastic' when ity corresponds to a real strain tensor 
|22j . In hyperelastic materials, on the other hand, the existence of quadratic free energy functional 

F = -^XtjkiUijUki (3) 
gives an additional symmetry under exchange of the first and second pair of indices, i.e. 

Xijki = Xkiij- (4) 

In the "pseudo-elasticity" theory, the vector Ui will be a novel geometric quantity (see below) , and 
the resulting tensor ity will be called a pseudo-clastic strain tensor. This tensor is still symmetric, 
as explained in section but the above additional symmetry is in general not present. 

We wish to construct general solutions of the equilibrium equations 

dio-ij = 0. (5) 

In order to close the problem for the stress tensor, a supplementary condition is needed which is 
the condition of compatibility, 

d 2 z u xx + d 2 x u zz - 2d x d z u xz = 0, (6) 

resulting simply from the fact that the tensor ity is built with the derivatives of a vector Uj. This 
relation does not depend on a specific interpretation of the tensor in terms of real deformations. 

The entries of the stress and strain tensors can be arranged in vector form, i.e. S = 
(o~ xx ,o~ zz ,<j xz ) T and U = (u xx ,u zz ,u xz ) T , giving a matrix representation of the elastic modulus 
tensor, 



S = AU, (7) 



where 





/ ^XXXX 


^xxzz 


2A, 


A = 


I ^zzxx 


^zzzz 


2A 




\ ^xzxx 


^xzzz 


2A. 



zxz 



(8) 



The factors of 2 are due to the symmetry under exchange of the last 2 indices of Xijki and uu- 
Now, we want to express the compatibility relation in terms of the stress tensor, so we need to 
express U in terms of S, i.e. 

U = 6S, (9) 
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where B = (Bij ) = A 1 . Then equation (jfJJ for an anisotropic medium is rewritten as follows 

n h ()-y J + BsySgEj - ■2ii h <),i).y., = o. (io) 

For an isotropic medium, Bu = i?22 7 B21 = -B12, B 3 i = B^ = 0, for i = 1,2, thus the equation 
reduces to A(a xx + a zz ) = 0. 

In the following, we will look for solutions of the form oc e lqx+luiz . j n this case, equation 
(|10f) . together with the conditions of mechanical equilibrium ijjjj, can be rewritten in matrix form : 

A{q,u>)X = 0. (11) 

A non-trivial solution occurs if det(_4(<7, lo)) = 0, which leads to a certain dispersion relation of the 
form tu(q) = Xq where X obeys the following equation: 

B22 _ B23 + 2B 32 ^ + 2B 33 + B 2 i + B u x2 _ gig + 2B 3 i ^3 + ^4 = q / 12 ) 
Bn Bu Bu Bu 

Depending on whether the roots X are real or complex, the response function will be qualitatively 
different: 

• Complex roots, corresponding to elliptic equations for the stress, appear within the classical 
theory of anisotropic elasticity. The fact that the roots are complex follows from the positivity 
of the free energy [23j . 

• Purely real roots can occur in the context of directed force chain networks considered below. 
The existence of at least one purely real root of the dispersion relation classifies the problem 
at hand as 'hyperbolic' [23I ]. 



B. The case of uniaxial symmetry 

Let us consider the case of uniaxial anisotropy and choose x and z to be along the principal 
axes of anisotropy. Then only Xijki with even numbers of equal indices is nonzero. Due to the 
symmetry of Xiju , this leaves one in general with 5 different constants. The matrix A takes 
the form 

a c 

A t = I J b ) . (13) 
d 

We denote it with a dagger to indicate that it corresponds to a material with a vertical uniaxial 
symmetry. An alternative parametrization of Aj, standard in elasticity theory, is 



1 - v x v z 




V * E * 



where E x>z and G are the Young and shear moduli respectively, and v X}Z the Poisson ratios. Note 
that the present form includes a linear elasticity theory without a free energy functional. The 
classical theory is recovered with the extra symmetry c' = c. In this case, E x , E z , v x and v z are 
not independent, satisfying the relation jjf*- = Together with G, we arc thus left with four 
independent constants. 

In classical elasticity theory for a uniaxial system, the stress-strain relation is derivable from an 
energy density of the form 

F = 5 [ au L + bu lz + 2cu xx u zz + 2du 2 xz ] , (15) 

The material described is stable under deformations if and only if F is positive definite for any 
strain, which requires 



a > 0; b > 0; d > 0; and ab - c 2 > 0. 



(16) 
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Or, equivalcntly, 

v x v z < 1; E x > 0; E z > 0; and G > 0. (17) 

An elastic material that is permitted to reversibly deform must obey these constraints, but they do 
not apply to materials for which there is no well-defined free energy quadratic in the strains. We 
speak of such materials as being described by coefficients that lie outside the "classical stability" 
range. 

The compatibility condition JHJl expressed in terms of the stress tensor reads: 

bd 2 a xx - cd 2 z a zz - dd 2 x a xx + ad 2 x a zz - 2^^-d x d z a xz = (18) 

Combining this relation with the two equilibrium conditions of equation (0 , 

d z a zz + d x a xz = 0, (19) 
d z a xz + d x a xx = 0, (20) 

we obtain, for any one of the components of the stress tensor: 

{dt + tdi + 2rdldl)a ij =Q, (21) 

where the coefficients t and r are given by 

a E x 
1 = b = £? 

ab - cd - \d(c + d) _ 1 / 2 v z v x \ 
" bd - 2 Ex \G ~ T. ~ E-J ■ (22) 

Expanding the stresses in Fourier modes, it is easy to see that the solutions of the equations 
(|19I21|) are of the form 

r*+oo 



/too 
dqY j a k {q)e^ X ^ z : (23) 
-°° k 

/+oo 
^^a fe ( g )X fe e wx ^ z , (24) 
-°° k 

/+oo 
dqY,M<l)X 2 k e^ x ^, (25) 
-oo t. 



where C xx and C xz are constants. From equation (|21|l we see that the X k arc the roots of the 
following quartic equation 

X 4 + 2rX 2 +t = 0, (26) 
a special case of equation (|T2*|) . There are four solutions: 



X = ±\l -r± x/r 2 - t. (27) 

Hence the index k runs from 1 to 4. The four functions ak(q) and the constants C xx and C xz must 
be determined by the boundary conditions, as shown in section ITTll and Appendix B. 

We see that only two combinations, r and t, of the 5 elastic constants will determine the structure 
of the response function in anisotropic materials. 



C. Main results of this paper 



We show in figure ^ the various 'phases' in the r-t plane corresponding to different shapes of the 
response function, as obtained from the calculation presented in section ITTT1 below. 



FIG. 1: (r,t) phase diagram characterizing the qualitative nature of the stress profiles. The shaded region 
corresponds to hyperbolic and "mixed" equations for stresses whereas the unshaded region allows for elliptic 
equations. The hyperbolic region is bounded above by the line t = r 2 , separating it from the elliptic region. 
In the elliptic region, a double peak stress profile is found in the whole region r < 0. The solid and dashed 
straight lines are the trajectories for the triangular spring network studied in section II VI for horizontal 
and vertical orientation of one of the springs, respectively. The symbols correspond to the solutions of the 
anisotropic linear dfcn model for various values of the anisotropic scattering parameter p, see section W\ 

The line t = r 2 , for r < 0, separates the hyperbolic and the elliptic regions. For t > r 2 (region 
I), the above roots Xk are complex and we write: 

X 1 = -Xt = /3-ia, (28) 
X 2 = -X 3 =-j3- ia, (29) 

where a and (3 are positive real numbers. When t < r 2 , r > (region II), one the other hand, the 
roots Xk are purely imaginary and one has: 

X 1 = -X 4 = -ia 1 , (30) 
X 2 = -X 3 = -ia 2 , (31) 

where ol\ and a 2 are positive real numbers. 

Note that the isotropic limit corresponds to the point r = 1, t = 1. As we show in detail in 
section lim the elliptic region contains a subregion r < 0, t > r 2 , where the response function has 
a two peak structure with peak widths growing linearly with depth. As one approaches the line 
t = r 2 , the two peaks become narrower and narrower, finally becoming two delta-function peaks 
exactly on the transition line. Below the transition, there is a hyperbolic regime (region III in 
figure ^) where the response consists of four delta-function peaks. 

The parameter range t < 0, labeled "mixed" in figure ^ gives rise to a third type of behavior 
of the response function due to the fact that there are two real roots and two imaginary roots. 
It may only appear in the non-stable pseudo-elastic case, and gives superposition of a hyperbolic 
two delta peak response function and a single-peka classically elastic response function. For the 
particular model for the DFCN discussed below, the range t < does not occur. Hence this case 
is not pursued any further here. 

We discuss below some particular trajectories in the r-t plane (see sections HVI and IV|) . One cor- 
responds to simple, anisotropic, triangular networks of springs, that lead on large scales to classical 
anisotropic elasticity with parameters on the plain and dotted straight lines, corresponding to two 
orientations of the lattice (see figure fTOf) . Both trajectories meet at the point (1, 1) corresponding 
to an isotropic medium where all springs have the same stiffness. Moreover, both trajectories cross 
the region r < and thus allow for two peak response functions. Inclusion of three-body forces 
permits spring networks with (r,t) anywhere in region I or II (see section llVjl . 
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We have also computed r and t for the linear DFCN model, for a particular model for scattering 
where the degree of anisotropy is tuned in terms of a parameter p (see section lv|) . The results are 
shown as symbols, and appear to always lie in the elliptic region. As in the spring networks, for 
sufficiently anisotropic scattering, one enters the region r < where the response function has two 
peaks. 

In two classical papers [24|, Green et al. have treated the stress distribution inside plates with 
two directions of symmetry with right angles to each other . The solutions are parametrized, apart 
from boundary conditions, by a\,a2 (not to be confused with cti introduced above) which are 
related to the set r,t by (r + \/r 2 — t/t), (r — y/r 2 — t/t). The authors assume their parameters 
Oil, &2 to be always real and positive, based on empirical fits of elastic constants for timbers such 
as oak and spruce. This choice corresponds to region II in figure^ Consequently, the possibility of 
region I and III behavior, and particularly the appearance of a double peak response for a classically 
elastic material, is not discussed in [24|]. Moreover, their analysis considers the response in the 
case where the boundaries and the directions of symmetry arc cither parallel or perpendicular to 
each other, whereas the present discussion - see in particular section III.B - treats a more general 
case. The response functions for region II, as computed in the present work, could in principle be 
reconstructed from the results of Y2M. 



After having discussed the general framework of anisotropic elasticity and the particular example 
of two-dimensional systems with uniaxial symmetry, we now turn to the actual shape of the response 
function in such materials. We will calculate the response of an elastic or pseudo-elastic slab of 
infinite horizontal extent to a localized force applied at the top surface. We shall consider the 
case of a semi-infinite system with a force applied at a single point on its surface, for which 
complete analytical solutions can be obtained. More general situations (finite spatial extension of 
the overload, finite thickness of the slab with a rough or smooth bottom, ...) should be considered 
to obtain quantitative fits of experimental |ld IT^ and numerical data. Still, two angles are left 
free: the angle 8q that the applied force makes with the vertical, and the orientation angle r of the 
anisotropy with the vertical. 



We are interested in the response of a semi-infinite system to a localized force at its top surface 
z = 0. We suppose that this force is of amplitude Fq and makes an angle 6q with the vertical 
direction, as shown in figure |21 

The corresponding stresses at z = are then 



III. SHAPE OF THE RESPONSE FUNCTION 



A. Vertical anisotropy 



Fo cos 9q S(x), 
Fq sin 0o S(x). 



(32) 
(33) 








> 



x 




FIG. 2: Force at the top surface. 
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FIG. 3: Region I: Rescaled stress profiles for several directions 80 of the applied force and several values 
of r, with t — 2. In each panel, the thick solid line is for r — —1.3, the thick dashed line is for r = —0.7, 
and the thin solid line is for r — —0.2, and the thin dashed line is for r = 0.5. r > is the condition to 
have a single peaked profile for 9q — 0. 



To obtain the results described below, we make use of the identity 

= rdq{e iq * + e- iqx ), (34) 
27T Jo 

and impose the boundary condition by identifying the coefficients of e ±lqx in the equations I|32I33[) 
and ^23 (1 -1)25 (1 at z = 0. Note that a xx (z = 0) is not determined by the boundary conditions. 

When z — > +00, we expect all stresses to decay to zero. It turns out that this is a self-consistent 
condition as long as the system is energetically stable, but cannot be imposed in the unstable 
regime. The reader interested in a more detailed derivation of the following results can consult 
appendix B. 



Region I (elliptic): t > r 2 

Since we want all the stresses to vanish at large depth, the functions 01 and a-i in 1)23 (1 -1)25 (1 must 
be zero for q > 0, and 03 and 04 must vanish for q < 0. In addition, C xx , and C xz must all vanish. 
Furthermore, because the stresses are real quantities, a\{— q) = a^(q) and a^{— q) = a%{q)- 

The boundary conditions at z = then imply 

F 

a 3 = —— [{13 - ia) cos O - sin 9 ] , (35) 
47tp 
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F 

a 4 = - — ~[(f3 + ia) cos 9 + sin 6>o . (36) 
Anp 

Since the coefficients 03 and 0,4 are independent of q, the integrals in equations (|23fl - l|25[l are 
straightforwardly carried out, yielding 

F 4az 2 [z cos 8 (a 2 + (i 2 ) + xsm9 ] 
° zz = 2^ [(a 2 - [3 2 )z 2 +x 2 } 2 + [2af3z 2 } 2 ' ( ^ 

(38) 



z 



2 



0- ra = - cr zz . (39) 



The latter two results follow directly from the observation that the integrals in equations (|24|) and 
(|25|l can be expressed simply as convolutions of a zz {q) with the Fourier transforms of x/z and 
x 2 1 z 2 , respectively. In the limit — > (which corresponds to r 2 — t — ► 0) and a — ► 1, we recover 
the familiar isotropic formulas [2l|. 

Figure shows the response for four different choices of the parameter r and a fixed t, each 
being shown for three choices of 9q. Note that a zz has a more pronounced double-peak structure 
for increasingly negative r. For #0 = 0, the condition for having a double peak is d 2 a zz (x — 0) > 0, 
which occurs when a 2 < /3 2 , or equivalently r < 0. In terms of the Young and shear moduli and 
the Poisson ratios, this condition can be expressed as G > E x /v x = E z /v z . The positions of the 
peaks are then given by x = ±zy/ j3 2 — a 2 = ±z-\/jr}. From the curvature at the maximum, one 
can define a width w of these peaks which reads: 



af3 1 Vt^V 2 , x 

w = -!=—= z= ■ z. (40) 

V2 ^2 _ a 2 ^y/W\ 

Thus the peaks become sharper and sharper as one approaches the hyperbolic limit t = r 2 . 

A very important point is that the response profiles scale with the reduced variable x/z when 
multiplied by the height z. This means that, when the profile is double peaked, these two peaks get 
larger in the same way that they get away from each other. Such a response cannot therefore be 
seen as an 'hyperbolic-like' signature, for which the peak width compared to the distance between 
the peaks goes to zero at large depth. However, in the limit where t — > r 2 , the width of the peak 
vanishes a the response becomes truly hyperbolic. 



Region II (elliptic): t < r 2 , r > 

Again, we only keep the functions a\ and ai for q < 0, and 03 and 04 for q > 0. This time, 
the fact that stresses are real quantities requires a\{—q) = 04(g) and a^(— q) = a^(q). A similar 
analysis to the above yields 

F 2(ai + a 2 ) z 2 {a 1 ot2Z cos 9 + xsm9 ] 
azz ~ 27 [(aiz) 2 + x 2 ][{a 2 z) 2 + x 2 } ' ( ' 

X 

Oxz = -cr zz , (42) 

z 

a xx = er zz . (43) 

For oc\ = a.2 = 1 (again r 2 — t = 0), we recover the isotropic formula. In this case, however, 
when 80 = 0, a zz always presents a single peak, see figure 0] Depending on the values of a\ 
and «2, the profiles can be broader or narrower than the isotropic response, as has been observed 
experimentally on, respectively, dense and loose packings 0| . 

1. Region III (hyperbolic): t < r 2 , r < 



In this case, all the roots Xk are real, and the response function is the sum of four 5 peaks, 
at positions x = X^z. The appearance of four peaks is different from previous hyperbolic models 
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-4-2 2 
rescaled horizontal position x/z 



FIG. 4: Region II. Stress profile for different cases. The solid thick line is for t = 1 and r — 1 (isotropic 
case), the thick dashed line is for t = 1 and r = 2.125, and the solid thin line is for t — 2 and r = f .5. 



giving two peaks in which case the closure relation for the stresses is linear, whereas here 
the closure is achieved by a 4th order partial differential equation, Ea. i|21|) . The four peaks merge 
into two peaks exactly on the hyperbolic-elliptic boundary t — r 2 . The reason why previous 
hyperbolic models 0, work so well could be that granular system such as sandpiles are close to 
the hyperbolic-elliptic boundary (see also section IV. B for further remarks). Inside region III, the 
fact that all roots are real excludes the possibility to require stresses to vanish for large z. This 
leads to a situation where there are more constants of integration than boundary conditions. 

One may advance on the analytical form of response functions using physical arguments as 
follows. Let us first rewrite the equation for stresses as follows 

(9f-c^)(^ 2 -c 2 _^K=0 (44) 

where 

4 = -r ± \jr 2 -t (45) 

leading to c± > 0. The constants ±c± are just the four real roots Xk mentioned above. Instead of 
solving the equation above, we consider special solutions er^, <j~j of the following PDE, 

(d 2 -4^)4=0, (46) 

which automatically satisfy equation (1441) . Both equations can be solved for the boundary condi- 
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tions i|32|) and (|33|l . giving the solutions 



cos 0q — 



sin Wo 



(5(x + c±z) + 



cos 6*o + 



sm Wo 



(5(x — c±z) 



2 

— (— [c± cos # — sin 8 ] S(x + c±z) + [c± cos 9 + sin 9 ] S (x 



c ± z )) ■ 



(47) 
(48) 



a xx = {c±[c± cos 6q — sin 9q]S(x + c±z) + c±[c± cos 6q + sm 9q]S(x — c±z)) . (49) 

Before constructing a general solution from tr- , let us remark that there are in principle additional 
solutions &ij satisfying 



[01 



c\dl)aij 



=F 



(50) 



However, these solutions are not finite as they involve divergences arising from integrals such 
as dq cos(qu) / q 2 . Therefore, we conclude that a general solution of equation l|44|) may be 
constructed as 



a+a. 



a_er„ 



It should satisfy the boundary conditions (|3"2")l and which yield a relation 

a + + a_ = 1. 



(51) 



(52) 



The coefficients a+ and a_ = 1 — a+ are relative weights which indicate how the applied load is 
shared between the two sets of force chains characterized by c± . As there is no physical mechanism 
introduced a priori which prefers one set of force chains to the other, we are left with one free 
parameter, say o+, for the response function c^. The ambiguity on the value of a+ could be 
resolved by considering e.g. a microscopic model that leads to equation (|44J) . 

In figure the propagation of the applied force along the characteristics is shown. Note that 
the sign of a zz may change along a certain characteristic if cosWq — S1 "_^ Q < (see figure EJ(b)). 



B. Anisotropy at an angle 



We now, for completeness, generalize the results of the previous subsections to the case where 
the direction of the anisotropy makes an arbitrary angle r with the vertical. (The previous section 
corresponds to r = 0). This situation may be relevant for systems that are initially sheared as in 
the experiments of Geng et al. |27j . or prepared in a way which breaks the symmetry x <-+ —x. 
We restrict the discussion to regions I and II (the computation for region III can be carried out in 
a similar fashion). 

The equivalent of the relation J7J) involves now a matrix A T which is related to Af of equation 
13 by 



a t = e _1 A t e, 



where Q is the rotation matrix 



. 2 

sm r 



Q = 



sin 2 r 



sm t cos r 



-2 sin t cos t 
-2 sin t cost 



sm t cos r cos r 



. 2 
sm r 



(53) 



(54) 



The differential equation on the stress components that is deduced from the compatibility condi- 
tion and stress-strain relations is now much more complicated, but the corresponding roots of the 
fourth order polynomial that appear when looking at Fourier modes can still be calculated from 
the Xk solutions of l|26|l . They read 



tanr 



1 + Xk tan t 



1, 



(55) 



The same method as above - see also Appendix B - can then be applied to find the stress response 
functions for a localized overload at the top surface of the material. Note that the material 
properties are still determined by the Xk associated with Aj. In particular, whether the response 
is elliptic or hyperbolic cannot depend on r. In the following, regions I and II arc defined with 
respect to Xk as above. 
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FIG. 5: Region III. The stress profile is a sum of four delta functions. The characteristics x — ±c±z along 
which the applied load is propagated are shown. Parameters are r = —1.0, t — 0.75 giving c+ = 1.5 (solid 
lines) and c_ = 0.5 (dashed lines). The delta functions are indicated by cartoons, (a) 8q — 0, (b) 8o — 7r/4. 



Region I 

The Xk are of the form ±/3 ± ia, see (|28I29|) . The corresponding can be constructed with 
the following quantities 

= a(l + tan 2 r) 

(l + /3tanr) 2 + (atanT) 2 ' k ; 



= g(l - tan 2 r) + tan r(q 2 + (3 2 - 
(l + /?tanr) 2 + (atanr) 2 



= a(l + tan 2 r) 

(l-/3tanT) 2 + (atanr) 2 ' V 7 



B , = /3(l~tan 2 r)-tanr(cv 2 +/3 2 -l) 
(1 - /3tanr) 2 + (atanr) 2 



The same boundary conditions - see figure [21 lead to 
Fn 2z 2 



{x sin 9q{A + A') (60) 



2tt [{x + Bzf + (Az) 2 }[{x - B'z) 2 + (A'z) 2 
+z cos 9 a [AA' (A + A') + AB' 2 + A' B 2 } + [x cos 9 + z sin 9 ](A' B — AB') } 
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FIG. 6: Region I. Response profiles for different values of the anisotropy angle r, but with a fixed value 
for the orientation of the applied force: 8o = 0. The graph (a) is for t = 0.6 and r = —0.2, while (b) has 
been obtained for t = 0.6 and r = 0.2. Note that for the three smallest r > the response is stronger in 
the negative x region. 



a xz and a xx are related to a zz by the usual factors of x/z and (x/z) 2 respectively. 

Figures El and show the pressure response profile as different parameters arc varied. In figure 
El the applied force is kept vertical (6q = 0), and r is varied from to tt/4. Interestingly, the 
initially double peaked profile (figure EJ- (a)) is progressively deformed in such a way that the left 
peak gets more pronounced, until the remaining single peak moves to the right for r = tt/4. This 
behavior might be counter-intuitive for smaller r, because a positive value of r means that the 
main direction of the anisotropy is oriented to the right. However, it can be understood within 
the ball-and-spring model of section Hvl where the hi springs are horizontal. Rotating to the right 
the two stiff directions hi emerging from a ball downwards brings the left one closer to the vertical 
direction, which therefore gets a larger fraction of the overload. Continuing past r = 7r/6, however, 
the stiffer springs form lines that slope downward to the right. Since they continue to support most 
of the load, the single peak is shifted to the right. This behavior holds also for the single peaked 
profiles of figure El(b). 

The second series of plots - figure [7]- is for the case where the applied force is exactly in the 
direction of the anisotropy (9q = r). The corresponding curves are qualitatively similar to those 
of Figure El The direction of the force imposed at the top does not change the general shape 
(anisotropic double or single peak) except for the fact that a negative pressure zone evolves for 
large negative x. 

The value of 0.6 for t used in the figures El and [3 is motivated by experimental findings |2q. The 
response function shown in figure E} (b) for r = 7r/4 is at least qualitatively consistent with the 
response functions measured in 27]. 



Region II 

In region II, where X\ = —X4 = —ia.\ and X2 = —X3 = —iot-2-, the expressions of the corre- 
sponding Yfc involve the quantities 

0^(1 + tan 2 r) 

Al = T~( 1 \2> ( 61 ) 

1 + (ai tanr) 2 

1 + {cti tanr) 2 
q 2 (l + tan 2 r) 

A ' 2 ~ T~ 1 — I \2> ( b6 

1 + (a.2 tanr)^ 

taiiT^-l) . , 

B 2 - -7—, 7 v?, (04) 

1 + (a 2 tanr) 2 
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FIG. 7: Same graphs as in figure |S] but this time with 6 = r as indicated in legends. 
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FIG. 8: Region II. Response profiles for different values of the anisotropy angle r, but with a fixed value 
for the orientation of the applied force: f?o = 0. The graph (a) is now for t = 2 and r = 1.5, while (b) has 
been obtained for t = 0.6 and r = 0.8. This time, the response peak can be moved to the right or to the 
left with positive values of r. 



the pressure response having the form 

Ft 2z 2 

° zz = ^Tu — , o \ 2 i (a — ^TTT — i p \2 i i a — {x sin 9 (A 1 +A 2 ) (65) 
27r [(x + B 1 z)' s + {A 1 zy\[(x + B 2 z) z + (A 2 z) z \ 

+zcos8 {A 1 A 2 (A 1 +A 2 ) + A X B\ + A 2 Bf] + [xcos6 Q + z sin 8 }(A 2 B 1 + A^)} . 
Again, the expressions of a xz and a xx are not shown, but can be deduced as usual from that of 

0~zz- 

The Figures Eland El show the response profile for different values of the parameters. Depending 
on these parameters, the response peak can be moved to the right or to the left with positive values 

Of T. 

Please note that the response function shown figure IHl-(b) for t = 7r/4 also agrees qualita- 
tively with the experimental findings in |27| . A more detailed analysis of their results is certainly 
worthwhile, also in order to possibly decide whether region I or II behavior applies for a sheared 
two-dimensional layer where the angle of the preferred orientation of force chains coincides with 

T = 7r/4. 
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FIG. 9: Same graphs as in figure but with 9 = r as indicated in legends. 




FIG. 10: Network of springs of stiffness fci and &2- 



IV. TRIANGULAR SPRING NETWORKS AND ANISOTROPIC ELASTICITY 

A. Triangular spring networks 

To illustrate the previous calculations, it may be useful to construct a ball-and-spring model 
with a tunable parameter that allows us to obtain different relative values of a, b, c, and d above. 
Here we consider a triangular lattice of balls with springs connecting all nearest-neighbor pairs. 
The lattice may be oriented in either of the two ways shown in figure 1101 and the springs have 
stiffnesses ki or k 2 as shown for the two cases. All springs lying along a given direction have the 
same stiffness. We take the equilibrium lengths of all springs to be unity. 

In either orientation, the system has reflection symmetry under x — ► — x and z — > — z, but not 
under rotations; it is described by an anisotropic stress-strain relation of the form of Aj . We deter- 
mine the elastic coefficients by writing down the energy directly for a homogeneous deformation. 
Note that the balls form a Bravais lattice, and hence that their displacements for a given average 
strain Uij are simply given by UijTj, where r is the equilibrium position of the ball. The energy 
density can easily be obtained by summing the energies of the three springs linking the ball at 
(0, 0) to its neighbors along different lattice directions and dividing by the area of the unit cell, 
A = V3/2. 

Horizontal orientation of the ki -springs 
For the case where the fci-spring is horizontal, we find for the energy density: 

F = — j— • [(8fci + k 2 )u 2 xx + 9k 2 u 2 zz + 6k 2 u xx u zz + 3k 2 {u xz + u zx f] , (66) 
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which corresponds to a matrix Af with the following coefficients: 

8ki + k 2 
a = —8A~' 
, 9k 2 
b = 8A' 

3fc2 

C ~ 8A' 
6k 2 

d = Si' 



(67) 
(68) 
(69) 
(70) 



Without loss of generality, we rescale all stiffnesses by a factor 8A/k2 and let k\jk 2 be denoted k. 
The coefficients r and t of equation 126(1 are then given by 

. - i±£, (7i) 

r = — , (72) 

which gives r 2 — i = ^-k(k— 1). We may eliminate fc from these two equations to obtain a trajectory 
in (r, t) space: 

(«) 

shown as the plain line in figure ^ 

Thus, fe < 1 (weak horizontal springs) corresponds to region I above with (see equation H28B ): 



, 2 



= - (4k - 1 + V8k + lj , (74) 



2 = - (1 - 4k + V8/c + lj . (75) 

As mentioned above, the condition for a double-peaked a zz profile is r < 0. Hence the single- 
peaked shape of o zz (x) becomes double-peaked when k < 1/4, i.e. when the horizontal springs are 
substantially softer than the others. 

For k > 1, on the other hand, we are in region II with (see equation (JSUJl): 



a\ = - [4k - 1 + 4yjk{k- l)j , (76) 



i (4k - 1 - Vk(fc-l)) • (77) 
The ct zz profile is always a single peaked when the horizontal springs are stiffer than the others. 



Vertical orientation of the ki -springs 

For the case where the fci-spring is vertical, we get a matrix Af where the coefficients a and b 
have been swapped from the horizontal case, i.e. with the following coefficients: 

a = iJ , (78) 
3fc 2 
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FIG. 11: Variables associated with three-body bond-bending interaction. 



Again, we rescale the stiffnesses and let k = fci/fo, this time finding 

1 = ITS' (82) 

3(4fc - 1) 

T = ~YT8k~' (83) 

which gives r 2 — t = ^^^^-P ■ As before, k may be eliminated to obtain the trajectory in (r, t) 
space: 



t = -2r + 3, (84) 



now corresponding to the dotted line in figure ^ 
For k < 1, we are in region I with 

« 2 = TTsh' < 85 » 

Again, the single peaked shape of the a zz profile becomes double peaked when k < 1/4. 
For k > 1, we have r 2 — t > and we are in region II, with 

a? = TTsfc ( 4fc ~ 1 + 4 v / M^ T )) ( 87 ) 

3 



12 ~ l + 8fc 



4fc - 1 - 4 



Three-body (bond-bending) interactions 

For the spring networks discussed above, the Poisson ratios are not both adjustable simultane- 
ously For the horizontal orientation of k± springs, v x = c/a is always 1/3, while for the vertical 
orientation v z = c/b is always 1/3. In order to have a ball-and-spring model on a Bravais lattice in 
which all clastic parameters can be varied independently, it is necessary to introduce three-body 
interactions. A straightforward way of doing this is to assume an energy cost for bond angles that 
differ from 60°. 

For simplicity, we present an analysis only of the case where the triangular lattice is oriented 
so that the k\ springs are horizontal. Consider the triangle of balls and springs shown in figure 
im We define 6y as ZXYZ, measured in the strained configuration. For the case of uniaxial 
symmetry, the energy of the triangle is determined by two bond-bending stiffnesses n\ and K2- For 
case I we define 



E bb = (1/2) 













-5)'- 


\-K2 ( 


»-!) 











K2 



(89) 



with Ki assigned to the angle opposite the horizontal edge. As for equation we take the 

equilibrium lengths of the springs to be unity. 
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Writing expressions for the angles in terms of displacements of the balls from their equilibrium 
positions and summing over all triangles, including the upside-down ones (shown dashed in figure 
on a homogeneously strained lattice, we find a contribution to the total energy density of 



"66 



3 

8A 



[(2ki + K 2 ){u 2 xx + u 2 zz ) - 2(2«i + k 2 )u xx u zz + 12k 2 u 2 xz ] . (90) 



Adding this contribution to equation l|6tj|) gives a total energy density corresponding to a matrix 
Af with coefficients 

8fci + fc 2 + 6K 

a = ' (91) 

9k 2 + 6k 

b = ^4—' 

3k 2 — 6k 

C = 8A ' 
_ 6(k 2 + 6k 2 ) 
~ 8A ' 

where k = 2«i + k 2 . In terms of bulk and shear moduli and Poisson ratios, we obtain 



(92) 
(93) 
(94) 



9fcifc 2 + 6(fci + 2k 2 ) K 
z ~ (8k 1+ k 2 + 6 K )A ' 
3k 1 k 2 + 2(k 1 +2k 2 )K 
E * = (3fc 2 + 2k) A ' (96) 

G = (97) 
3fc2 6K (98) 



8fci + k 2 + 6k 
k 2 — 2k 
3fc 2 + 2k 



(99) 



Note that E x v z = E z v x , as expected. Note also that it is not necessary for k\, k 2 , K\, and k 2 to 
all be positive. Stability (cf. equation lfT7)l) requires only 

8/ci + k 2 + 6k > 0, (100) 

3fc 2 + 2 K > 0, (101) 

3fcifc 2 +2K{h +2k 2 ) > 0, (102) 

k 2 + 6k 2 > 0. (103) 



From equation ill'L'li we find 



3(3/e 2 + 2k) 

^ _ /4\ 3fcxfc 2 + 2 K (h + 2fc 2 ) 4fc 2 

(3fc 2 + 2 K )(fc 2 + 6K 2 ) 3fc 2 + 2ft' 1 j 

By choosing ki, k 2 , and k we can obtain any positive value for t. From Eqs. i|lUU|) - (|103f) . we see 
that the second term in the expression for r is positive. For fixed t, we can make r arbitrarily 
large by choosing k 2 close to —k 2 /6. The smallest (or largest negative) value of r is obtained by 
choosing 3fcifc 2 + 2n(kx + 2k 2 ) = (and adjusting fe 1; say, to keep t fixed). This leads to r 2 — t = 
and r < 0, demonstrating that the triangular lattice can lie anywhere in region I or II. 



B. Remarks 



We have seen that classical anisotropic clastic materials can have double-peaked response func- 
tions and that such cases can be obtained with simple ball-and-spring models. These calculations 
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explain, for example, the numerical results of Goldcnbcrg and Goldhirsch [l8j, without invoking 
any special considerations on small system sizes. 

It is important to note that the response functions for the triangular spring networks always lie 
in the elliptic regime: the peaks broaden linearly with depth. Thus the observation of a double- 
peak structure is not necessarily an indication of propagative (hyperbolic) response in an elastic 
material. However, when the k\ springs are oriented horizontally, and in the limit where their 
stiffness tends to zero, the response becomes hyperbolic. In this case, one generically expects 
peaks to broaden diffusively, i.e. like \J Dz 0, [2f|. Note that in the limit where ki — > 0, there 
appears a floppy (zero energy) extended deformation mode which, as emphasized by Tkachcnko 
and Witten [l(| , naturally leads to a stress-only closure equation and hypcrbolicity. In the phase 
diagram, figure ^ this limit corresponds to the point where the straight solid line touches the 
boundary curve t = r 2 . Note that within this line of thought, one should also expect hyperbolic 
response in elastic percolation networks at the rigidity threshold. In fact, in the limit k\ — ► the 
triangular network becomes a rhombic network which is known to become isostatic for a finite 
system: a single boundary suffices (say a bottom surface in the slab geometry) in order to suppress 
the zero mode, and the system becomes rigid @ . 



V. ANISOTROPIC DIRECTED FORCE CHAIN NETWORKS 



A. Biased scattering 

In |13| . a Boltzmann equation for the chain-splitting model was derived for a granular medium 
which is strongly disordered. In the present work, we suppose that the scattering of force chains by 
defects is biased by a preferred orientation of the material, modelled in terms of a global director N. 
We intend to describe systems possessing a uniaxial symmetry which have undergone compaction 
or shearing or which have been constructed by sequential avalanching due to grains poured from 
a horizontally moving orifice. 

The fundamental quantity is the distribution function P(f, n, r), where: 

P{f,n,r)dfdnd D r (106) 

gives the number of force chains with intensity between / and / + df, inside the (solid) angle dn 
around the direction n, in a small volume element d D r centered at r. Integration of P(f, n, r) with 
respect to / and n will yield the density of force chains at the point r. [3(j The distribution function 
is defined with respect to an ensemble of different realizations of force chains for an assumed 
uniform spatial distribution ofpoint defects (of density p^), with same boundary conditions. In 
the spirit of previous models [jj, |26j that give hyperbolic equations for the stresses, a mechanism 
of propagation is implemented, but now on the local level of force chains. In the analytical model 
presented here, a pairwise merger of force chain to a single one will be neglected. The limitation 
of this approximation will be discussed below. Then the distribution function P(f, n, r) obeys the 
following linear equation 

P(/i > ni,r + n 1 dr)= - P(/i,n 1)r ) 

+ 2y Jdf jdf 2 jdrx' /dnaP(/V,r)#(n'->ni,n a |N) 

x 5(/ 1 cos^ 1 +/ 2 cos^-/ / )5(/isine 1 + / 2 sin^)|sin(ei-e 2 )|, (107) 

where A is the mean free path of force chains, and is of the order of 1/ '(pdl D1 ) in D dimensions. 
The length I represents the average size of a grain. The equation means the following: a force 
chain at some point r + \i\dr is either due to an unscattered force chain, which occurs with the 
probability that no scattering occurs times the probability that the same force chain existed at 
point r (given by the first term on the r.h.s. of the equation), or to a scattered force chain. The 
latter occurs with the probability given by the second term on the r.h.s. of the equation: it is 
the sum with respect to all intensities and directions of the incoming (labeled by a prime) and 
the second outgoing force chains of the product of the probability for the incoming force chain to 
arrive at r times the probability of scattering ^^(n' — > ni,ri2|N). The delta functions impose 
conservation of forces, the factor 2 accounts for the number of outgoing force chains, and the 
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factor | sin(#i — 6%)] is convenient to write explicitly rather than include in ^> . The dependence of 
the scattering probability on N requires to consider the outgoing force chains separately. In the 
absence of N the outgoing force chains may be treated symmetrically, and one recovers the linear 
model for an isotropic medium [l3j . 

The analytical model presented for biased force chain scattering does not take into account fusion 
of force chains, which leads in general to a non linear Boltzmann equation. For an isotropic medium 
the consequences of fusion have been discussed for a model where force chains arc restricted to 
lie on exactly 6 directions In this discrete model the validity of the linear approximation 

was explicitly shown to be restricted to shallow systems (depths smaller than a few times A) and 
small forces. However, preliminary results on a discrete model with 8 directions suggest that the 
linear theory might have a wider scope of application than expected from the study on the 6-leg 
model. More precisely, a proper analysis of the linear perturbation analysis around the full non- 
linear solution of the Boltzmann equation might share, in some regimes, many properties of the 
linear solution presented here. In any case, one can see the present analysis as a shallow layer 
approximation where the fusion of chains can indeed be neglected. 

Instead of solving equation i|l(J7|) . we first introduce the scalar local average force density F(n, r), 
i.e. the local scalar force field per unit volume, defined as 

y>OC 

F(n,r) = / df /P(/,n,r). (108) 
Jo 

Then, multiplying equation I|1U7|) by /, we obtain the following equation for F(ni,r): 

Ani • V r F(ni,r) = -F(ni,r) 
+ 2 Jdv! fdn 2 F(n',r)^(n' — m,n 2 |N) 
1 

cos 9i — (sin 0i / sin 6*2 ) cos 8 2 ' 

This equation is identical in form to the Schwarzschild-Milne equation for radiative transfer [3l| . 
though, unlike the situation in radiative transfer problems, the albedo is larger than unity. Let 
us note that the possibility to rewrite the Boltzmann-type equation I|1U7|) in terms of the force 
density F(n, r) is only possible for the linear model. 

/From now on, we take all lengths in units of I which amounts to formally setting 1 = 1. Now, 
we introduce physically relevant angular averages 

p(r) = JdnF(n,r) (110) 

Ji(r) = JdnmF(n,r) (111) 

ofc(r) = D J dnriiTij F(n, r), (112) 

where J dfl is a normalized integral over the unit sphere. The field p is the isostatic pressure, 
while J may be interpreted as the local directed average force chain intensity per unit surface. 
Now, given a local snapshot of a force chain network, one can usually not tell the direction of each 
chain. Moreover the average force vanishes everywhere in the system as a consequence of Newton's 
third law. The directions of chains arc actually determined by the boundary conditions, say on 
the top and bottom of a granular layer, which thereby determine the field J in the bulk. It is the 
propagation of force chains starting from the boundaries of the system modeled by equation (|107fl 
which leads to the orientation of the force chain network. Finally, the tensor a is the stress tensor. 



B. Stress equilibrium at large length scales 

We now proceed to obtain the equations governing the physically relevant fields introduced 
above, by calculating the zeroth, first, and second moment with respect to rij of equation (|lUtJ|> . 
The equations read as 



A V- J = (ci - l)p + c 2 (J NN , 



(113) 
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djtra = 0, (114) 



(5y V • J + 0jJj + djJi) = B Q o-j 



(D + 2) 

+ S tJ (Bi AV • J + B 2 a NN ) + NiNj- (B 3 AV • J + B 4 a NN ) 

+ B b (N l( J 3k N k + N 3 <j lk N k ). (115) 

where (J n n = N • er • N. The second equation (|114() is readily obtained upon averaging, while the 
first and third, equations (jll3|) and ({115(1 , are obtained using an Chapman-Enskog-type expansion 
of the local average force density F(n, r) in terms of the fields p, J, and a already given in [lflj : 

D + 2 

F(n, r) = p(r) + Dn • J(r) H — n • o-(r) • n + . . . (116) 

Let us remark that equation 1114|) gives mechanical equilibrium as expected and is independent on 
the specific form of ^(n' — ► ni,n2|N). The validity of the Chapman-Enskog expansion is based 
on the assumption that on large enough length scale an isotropic state is reached. For the case of 
biased scattering of force chains considered here, this implies that the bias intensity must not be 
too strong. Then the statistical weight of the set of force chains propagating through the entire 
system without changing their direction will not be important. The limiting case of strong bias 
requires a different approach than the one presented here. 

The constants and B^ appearing in Eqs. I|113[) and l|115|) respectively are angular integrals 
involving the microscopic model for the scattering process, i.e. a specification of ^(n' — > n.i, n 2 |N). 
A specific model will be considered in the next section. If one neglects the dependence on N in 
the equations ab ove, one recovers the simpler equations for force chain splitting in an isotropic 
granular medium |l3j . 



C. A linear pseudo-elastic theory 

As in the isotropic case, one would like to see if equation (|115f) can be cast into a form where 
the stress tensor cry is a linear function of a pseudo-strain tensor 



djJi) , 



giving rise to the relation 



= A 



ijklUkl, 



(117) 



(118) 



where Xij k i is the anisotropic pseudo-elastic modulus tensor. Similarly to conventional elasticity 
theory as mentioned in section [HI we will see that the tensor Xijki satisfies the symmetries given 
in equation (J2J) . The symmetric form of u%j stems from the symmetries appearing in the derivation 
of the large scale equations when carrying out angular averages, in particular J dnniUjUkni. In 
equation (|115fl . the gradients of the field Ji appear only in combinations such as V- J and diJj+djJi. 
Please note however that unlike in classical anisotropic linear elasticity theory, in the present case, 



Kjkl 7^ ^ 



klij 3 



(119) 



except for certain cases imposed by the details of the scattering process. The absence of the 
symmetry present in the classical theory is possible because there is no underlying free energy 
functional. 

The relation between the stress tensor and the pseudo-clastic strain tensor can be derived using 
the second moment equation (|115|l . The latter can be rewritten in the following form 



J a — Bij k i<j k i, 



(120) 



where 



AV • J 



1 

D + 2 



Bi 



BzNiN.; 



A 

D + 2 



{diJj + djJi), 



(121) 
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and 



Bijki = -^-{SikSji + Siidjk) + B 2 S lj N k Ni, 



B5 
2 

BiNiNjNkNi. (122) 



+ -^(fyi\yv fe + S jk NiNi + SuNjNi + SuNjNk) 



The relation between J^- and Uki can be inverted to give 



a ij ~ 7;AijkiJki, (123) 



where Aijki has the same form as Bijki with the constants -B M being replaced by constants 
which arc obtained from the relation 

In particular, one obtains the following relations for the constants A^: 

Ao = (125) 
Bq 

2_Z?2 

M = ~ B (B + B 2 + B 4 + 2B 5 ) 7 (126) 

4 5 (Bo + S2 + 54 + 2B 5 ) ' 1 ' 

Bo(Bo + -B5) ^ ^ 

Now, one can finally determine the pseudo-clastic modulus tensor in terms of the tensor A^V- 

Kjki = — ^ f^yM + — — (BiAij mm + B 3 Aij mn N m N n ) Ski- (129) 

Thus, the pseudo-clastic modulus tensor Xijki becomes - via the tensor A^f-i and the constants A^ 
- a function of the constants B^ which depend on the specific scattering model used. 

In the next section, a special case will studied which allows us to derive a simple, but non-trivial 
equation for the stresses which supplemented by the mechanical equilibrium condition (jll4|) opens 
a way to determine the stress tensor, or, put differently, the response function. 



D. A microscopic model for force chain splitting in presence of a bias 

As mentioned in the previous section, the entries of the pseudo-elastic modulus tensor depend 
on the specific model for anisotropic scattering which is specified in terms of the scattering cross 
section conditional on the global director N, ^(n' — ► ni,n2|N). We have considered a specific 
model for force chain splitting. It tunes the strength of the bias for scattering parallel to N, using 
a weight for each outgoing chain proportional to powers of a cosine factor quantifying the degree 
of collincarity with the global director N (sec figure IT2"|) . 

For each force chain arriving at a defect in the direction n' two outgoing force chains are chosen 
in the directions ni and 112 as follows: the angle of one chain, say number 1, w.r.t. the incoming 
force chain is chosen with weight oc (ni ■ N) 2p , for a positive integer p, in the interval [0, 6 m ax] (or 
[-S mall 0]) while the other outgoing chain, say 2, is chosen uniformly in the interval [—9 max , 9\\ (or 
[— 0i, 9 ma x] respectively). The reason for choosing the direction of the second chain like this is that 
the first (biased) chain should carry most of the intensity of the incoming force. Increasing p leads 
to scattering which is more and more biased in the the direction N. The form of the scattering 
cross section is therefore chosen as 



*(n' m,n a |N) = C p ty(0 a |0i)(ni ■ N) 2 ? + ^(0i|0 2 )(n 2 ■ N) 2p ) 



(130) 
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N 



-e 



max 



FIG. 12: The microscopic scattering model. The length of the arrows are different to illustrate the amount 
of force transmitted along the directions. 



The functions ijj(6i\6j) are the respective (uniform) probabilities for 9i given Oj described above. 
The constant C p is a normalization factor which depends on the angle between n' and N and which 
is determined from 



and its explicit form is given in the Appendix A. 4. 

The simplest choice for the global director is N = z, i.e. if force chains are scattered preferably 
downward. We might think of a granular layer that has undergone compaction by a vertical load. 
In this case, the matrix Af relating the stress and pseudo-strain tensor has the block-diagonal form 
as given in equation (|13|) . Any other orientation of N can be related to the vertical one by an 
appropriate rotation (see section ITlIBI equation 153(1 ). 

The numerical values of the parameters r and t that determine the shape of the response function 
(see figure ^) depend in the case of the anisotropic linear directed force chain network model on 
the constants introduced in the previous section. The latter are calculated from the above 
microscopic scattering model (see Appendix A) and are listed in Tab. I-IV of Appendix A for 
different choices of the maximum angle max of the scattering cone and different bias intensities p. 

Interestingly, the roots we find for this scattering model all lie in the (elliptic) regions I and II 
introduced in figure ^ Hence, it is possible to find an anisotropic scattering rule that leads to a 
two-peak structure of the response function, but in no cases the values of r and t have been found 
to lie in the hyperbolic region. Whether this is a limitation of the linear treatment of the DFCN, as 
suggested by the analysis of the 6-fold model is at present not settled. Work in this direction 
is underway |32T |. 

We finish this section with the following remark. If one identifies the elastic constants of classical 
anisotropic elasticity theory and their geometrical generalizations obtained for the linear anisotropic 
DFCN, as we have always done implicitly here, the possible range of values which occur for typical 
granular materials can be discussed. Experiments indicate that in samples of sand which are 
filled from above and where the major principal axis of a stress tensor is in the vertical direction, 
t = E x /E z attains values in the range 0.4 < t < 1 (see H3]). For the maximum scattering angles 
plotted in figure ^ the values of t determined from the specific microscopic model for biased force 
chain scattering used here appear to satisfy the experimental range. Further information on the 
construction history of the sand samples, which affects e.g. the distribution of packing defects or 
the strength of the scattering bias, is needed to fully judge the quality of the anisotropic DFCN 
model presented here. 




(131) 



VI. CONCLUSION 



The main objective of this paper was to work out in details the response function to a localized 
overload in the case of linear anisotropic elastic, or pseudo-elastic materials in two dimensions. 
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After working out the details of two specific microscopic models, a triangular network of springs 
and an anisotropic directed force network, we have shown that the resulting large scale equations 
can lead to a large variety of response profiles, summarized in the phase diagram shown in figure 
n spanned by a two-parameter combination of entries of the (pseudo-)elastic modulus tensor. 
The one peak structure of conventional (elliptic) isotropic elasticity can split into two peaks for 
sufficiently anisotropic materials. This situation occurs as soon as the shear modulus G is greater 
than the ratio E x /v x — E z /v z of the Young modulus and the Poisson ratio (either in vertical or 
horizontal direction). This corresponds to an anisotropic material for which vertical stresses are 
easily transformed into horizontal strain (large Poisson ratios) and vice versa but which strongly 
resists shear stresses. However, contrarily to the prediction of 'stress-only' hyperbolic models, these 
two peaks generically spread proportionally to the height of the layer, and not as the square root 
of the height for an hyperbolic medium. For the triangular network of springs, there is a special 
point, where the lattice loses its rigidity and a soft mode appears, where the system becomes 
exactly hyperbolic. It would be interesting to exhibit other situations where these extended soft 
modes discussed in |lOj naturally appear; a possible candidate is a percolating network of springs 
at rigidity percolation. 

For the anisotropic rules of force chain scattering that we have chosen, on the other hand, the 
directed force network was always found to be in the elliptic regime. This might however be 
an artifact of the linear approximation that we have used and where mergers of force chains are 
ignored. Preliminary results suggest that for the full non-linear problem, a genuine elliptic to 
hyperbolic phase transition might take place when the degree of anisotropy is increased, but more 
work (underway) is needed to confirm this potentially interesting result. 

Recent experiments |27j have not been able so far to distinguish between a noisy hyperbolic 
response (where the width of peaks scales as the square root of the height) or anisotropic (pseudo-) 
elastic response functions. For sheared system where force chains are preferably oriented at 45 
degrees with respect to the vertical, response functions show a horizontal shift (in the lateral 
direction with respect to the point of applied force) of the maximum, consistent with the preferred 
orientation of force chains. We found qualitative agreement with our findings. More detailed 
experiments appear to be necessary to decide on the parameters r, t, i.e. the possible locations 
in the phase diagram, figure ^ or P U1, differently on the elastic constants, corresponding to a 
particular form of the response function, if the present (pseudo-)elastic analysis applies. 

It would be interesting to extend the present results to three dimensional situations in order to fit 
the results of experiments on deep sand layers, where a single peak response function was measured 
fl6|| . and, most importantly, to test the consistency of the effective elastic moduli obtained from 
this fit in other geometries (like the sandpile or the silo). It would also be very interesting to find a 
way to prepare a disordered granular medium in a sufficiently anisotropic state such as to observe 
a two-peak response functions. 
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Appendix A: Some integrals for biased linear dfcn 
A.l. Zeroth moment 

First, we propose to calculate the coefficients ci and C2. Using the expansion (|1 1 f>|> the integral 
w.r.t. ni of the equation for the force density, one finds 



AV-J = -p + 2 J dri J dn, J dn 
x *(n' -> m,n 2 |N) 



D + 2 

p + D ri i J l H — n^&ijrij 



cos 9 1 — (sin 6 1 / sin 62) cos 82 



= (k!-l)p + k 3 — ^a NN , (132) 

Please note that a contribution occurs only from terms which are even w.r.t. n' — > — ri . The first 
coefficient is given by 

Id = 2 [dri [ dn, [ dn 2 *(n' ^ ni ,n 2 |N) —J—— (133) 

J J J cos 0i — (sin 0i /sin #2) cos 02 

The second coefficient ^3 appears when performing a decomposition of the tensor 

2 [da! [ dn, [ dnanJn^n'-^m.nalN) . 1 . - 

J J J J cos 0i — (sin 0i / sin 02 J cos 02 

= fc°<5y + fc 3 7VjiV 3 - (134) 

The coefficient fcg is irrelevant because Sij&ij = where o-y = cr.y — <5yp. The the coefficient /C3 is 
given by 

fc 3 = 2 J dn' J dn, J dn 2 (2(n' ■ N) 2 - l) *(n' -> ni,n 2 |N) 

(135) 



cos 0i — (sill 0i / sin 2 ) cos 2 
One finally obtains 

D + 2 D + 2 

ci=h — fc 3 , c 2 = — - — k 3 (136) 

Explicit expressions for the constants ki, k 3 are given in section B.4 which finally will have to be 
evaluated numerically. 



A. 2. First moment 

Next, let us derive the equation of mechanical equilibrium pi4fl . Taking the first moment of the 
force density equation without an external force gives 



D 



= + 1 j dri J dn, J dn2n lti F{ri, r)*(n' -> m, n 2 |N) 



1 



cos 0i — (sin 0i / sin 2 ) cos 2 
The second term contains the integral 



(137) 



/ dn, [ dn 2 n hl -$(ri -> rn,n 2 |N) 1 — = ari t (138) 

J J cos 0i — (sin 0i/ sin 02) cos 02 

Symmetrizing the integrand w.r.t. to the indices 1 and 2 gives a = 1/2. This result is independent 
on the specific form for the scattering cross section ^(n' — * ni,n.2|N). The remaining integral 
w.r.t. n' yields J; canceling the first term — Ji above. 
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A. 3. Second moment 



Finally, we calculate the coefficients in the third of the hydrodynamic equations, equation 
(|115f) . Let us consider the second moment by multiplying the force density equation by ni^nij 
and integrating w.r.t. ni. One obtains the following equation: 

XDT m d k Ji = -lofc + J dn'F(n', r)I l3 (n', N) (139) 



where 



I y (n',N) = 2 J dm J daznitn^ip! -► n lf n2|N) 



cos 8i — (sin 6i / sin 62 ) cos 62 



(140) 



and 



^ijkl = + 2) ^ i ^ kl + ^**^3'i + ^jfc) (141) 

The tensor ly may be decomposed as follows: 

Iii (n', N) = A % + Kw'in'j + K 2 (n^Nj + n^-JV*) (142) 

The coefficients Ao, -Ki, an d ^2 are all functions of the argument n' • N which will be suppressed 
in the following. As the tensor Lj(n',N) should be invariant w.r.t. to the operation N — ► — N 
because the scattering cross section ^(n' — ► ni,n 2 |N) is, the functions K and K x are even and 
K 2 is odd under this "parity" change. They are to be determined by multiplying ly as follows: 

Iq=Iu = K D + K 1 +2K 2 (n' -N) (143) 
I x = NilijNj = K + Ki(n' ■ N) 2 + 2K 2 (n' ■ N) (144) 
J 2 =<T^ = A' + K\ + 2A 2 (n' ■ N) (145) 

(146) 

The variables Iq, I\, and I2 are likewise functions of the argument n' • N which is suppressed 
henceforth. In the following we consider D = 2. The system of equations may then be written in 
matrix form 

{h,hJ 2 ) T = A(A ,A!,A 2 ) T (147) 

with 

2 12 cos a 
A = ( 1 cos 2 a 2cosa ] (148) 
112 cos a 

where cos a = (n' • N). We eventually want the functions A M as a function of the integrals 1^ 
(i — 0, 1, 2. So we need the inverse matrix 

10 -1 

l/sin 2 a l/sin 2 a (149) 

-1/(2 cos a) 1/(2 cos a sin 2 a) — cos(2a)/(2 cos a sin 2 a) 



We find 



K = h-h (150) 

Ki = -r^-(I 2 -h) (151) 
sin a 

K 2 = — !— (-Iq + -^2—{h - cos(2a)/ 2 ) ) (152) 
2 cos a V sira / 
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Before writing down the integrals / M , let us introduce the vector 

N - (n ■ N)n 

n± ~V^ ( n ' N ) 2 
Furthermore, let us denote the integrals 



(153) 




(ni • n') 2 

(a) = 2 /dm /dn 2 ( (m-n^) 2 ] *(n' m, n 2 |N) 

J J * (m • 0(m • n') / 



(154) 



cos 6*i — (sin 6 1 / sin# 2 )cos# 2 



Then, we find for the functions A M 



K = h (155) 
K\ = —i\ + io — 2 sgn(sin a) cot m 2 (156) 

A' 2 = sgn(sina) (157) 

sin a 

The transformation from I , J 1; I 2 to io, i±, i 2 is primarily for technical reasons as in the scattering 
function the directions ni and n 2 are parametrized w.r.t. n'. We may now proceed to perform the 
integral on the r.h.s. of equation 11391) 

j du'Ffn', v)I t3 (n', N) = j da' (p + DJ k n' k + ^±la u n' k n^j Iy (n', N) 
= pO>H Jdn'K + j ' dn'K 1 n' i n' j + j dri 'K 2 {n' i N j + n^N,)^ 

J dn' Kin' i n'jn' k n' l 



D + 2 



&kl [ Sij I d-n!K rik n 'i 



2 

+ J dn'K^Nj + n'jNiWtri, 

(158) 

The integrals which are multiplied by J k give no contribution because due to their tensorial prop- 
erties they should all be linear in N which means that they are uneven under sign change N. On 
the other hand, the integrands are even w.r.t to this operation, which implies that the integrals 
are zero. 

We now further simpflify the integrals w.r.t. n' multiplied by p and a k i using decomposition 
according to Cartesian tensors. The integrals following p are denoted as follows 



dn'K = K (159) 

J dn'Km'rt = k lta Sij + K ltb NiNj (160) 

dn'K^Nj + n'jNi) = K 2 J VJ + A^A^- (161) 
The constants are given by 

K ha = / do;Aisin 2 a (162) 



Ki, b = J daKi cos(2a) (163) 
A 2 , a = (164) 
A 2 , 6 = 2 / daK 2 cos a (165) 
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The angular integrations above (and all the ones following below) are understood to be normalized 
by factors \/(2ir). The integrals following a 1 $ are the following: 

J dn'Kon'^ = + K^NiNj (166) 

M ijk i = J dri Kin' i n' :j n' k n' l = K 1 (<%4i + kk&ji + 5a8 jk ) 

+ K 2 (NiNjSki + perm. . . .) 

+ K 3 NiNjNkNi (167) 

and 

J dn'K 2 (n' i N j + n' j N i )n' k n' l 

= K' 2 (2N i N j S kl + NiNkSjt + N 3 N k 8 a + NiN,S jk + N 3 NA k ) 

+ 2K%N i N j N k Ni (168) 

Let us turn to the first of these three integrals. The coefficients are given by the following integrals 

K$, a = / daA'osin 2 a (169) 



K ,b = J daK cos(2a) (170) 

(171) 

The second integral l|167f) giving rise to the coefficients is treated by performing the following 
contractions: 

Mi = Mujj = J dn'Kt (172) 

= M '- M = h Kl cos2 ° (173) 

M 3 = M ijkl N i N j N k N l = J dances 4 a (174) 

(175) 

In matrix notation, the system of equations we have to invert is the following: 

(M 1 ,M 2 ,M 3 ) T = B(K u k 2 , K 3 f (176) 

with 

/ D(D + 2) 2D + 4 1 \ 
B = D + 2 D + 5 1 (177) 
V 3 6 lj 

Then, for D = 2 one finally obtains for the coefficients 

Ki = J daKi Q - ^ cos 2 a + i cos 4 a^J (178) 

K 2 = I daKi ( — — + — cos 2 a — — cos 4 a j (179) 



K 3 = J doK x (1-8 cos 2 a + 8 cos 4 a) (180) 
Finally the coefficients of the third integral l|168|) read as 

K' 2 = J daK 2 cos a sin 2 a (181) 

K 2 — [ daK 2 cos a(l- 4 sin 2 a) (182) 
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Next, one collects all coefficients in front of the Cartesian tensors on the r.h.s. of the second 
moment (113911 . 



XDT. ljkl V k Ji = --a VJ + aij{D + 2)K 1 



1 

D 

+ Sij (a p + aia NN ) + N t Nj (a 2 p + a 3 a N N) 

+ a A (Nja ik Nk + N t a jk N k ) (183) 



where the coefficients a„ are given as follows: 



ao 


= K + K ha 




(184) 


Ol 


D + 2 (k 


+ K2) 


(185) 


"2 


= K lib + K 2ib 




(186) 


«:i 


- D + V3 + 


2K<<) 


(187) 


0,4 


- (D + 2)(K 2 


+ K' 2 ) 


(188) 



and 



where 



(189) 



When reducing the integrals in terms of the integrals i^, one obtains 

ao = J doe (io sin 2 a + i\ cos 2 a — i 2 sgn(sin a) sin(2a)) (190) 
ai = — ^— dai 2 sgn(sino!) cos(2a) + K 2 ^j (191) 



a 2 = J da ((io — i\) cos(2a) + 2i 2 sgn(sina) sin(2a)) (192) 
D + 2 f 

a 3 = / da ((i - h)(l - 2sin 2 (2a)) 



2 

+ 4i 2 sgn(sina) sin(2o!) cos(2a)) (193) 

a 4 = (D + 2)(K 2 + i J da i 2 sgn(sin a) sin(2a)) (194) 

(195) 

K\ = J da ((io — ii) sin a — 2i 2 sgn(sino;) cos a) (_ 1 _|_ 4 cos 2 a) (196) 
~ / sin Ol 

K 2 = / da ((i — i\) sina — 2i 2 sgn(sina) cosa) — — (— 1 + 4cos 2 a) (197) 

Using the equation for the zeroth moment 1113(1 , we can eliminate p and we obtain the coefficients 
Bo through B$ 

Bo = -^ + (D + 2)K 1 (198) 
Bi = . 1 -r (a - 01 - (D + 2)^i) (199) 



B 2 = 01 + — ^— ( -ao + 01 + (D + 2)JTi ) (200) 



(ci - 1 

~i- 1 

B 3 = -^-r (a 2 - a 3 - 2a 4 ) (201) 
ci - 1 

Bi = a 3 + (-a 2 + a 3 + 2a 4 ) (202) 
ci - 1 

B 5 = a 4 (203) 
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Inserting for Ci, c 2 calculated in section B.l, and for K±, K2, and ao through 0,4 given above, the 
coefficients are entirely determined in terms of integrals over the scattering function or in terms 
of integrals over the functions z M which have to be evaluated numerically. Explicit expressions for 
the functions for a specific scattering model are given in section B.4. and have been used to 
yield the following Tab.s (I-IV) in the main text. 



A. 4. The scattering model 

We give now the explicit form for the normalization factor C p of the microscopic scattering 
model, equation l|13U|) . Choosing the angle between n' and N as a one finds the following relation 
to determine C„: 



I = dn x 



Cp2 



dn 2 *(n' -> m,n 2 |N) 

- do 1 r 01 do 2 



J max J — 



(cos 2p (0i 



a) + cos 2p (6>i + a)) 



(204) 



One finds 



c P {a) = 3 



J_(2p 

2 2p I p 



p-1 

E 

k=0 



2p \ sin((2p - 2k)9 max ) cos((2p - 2k)a) 



-, -1 



k 



(2(P-*)) 



(205) 



We have mentioned above that all constants of the hydrodynamic equations depend on the pa- 
rameters ki, and the integrals of the functions i M (a) for fi = 0, 1,2. Using the model for the 
scattering cross section introduced in the main text, they read as follows 



= 2 



„ (2^)° piaj l v cos(2a) / 



d6o 



,2p/ 



a) + cos 2 P(0 1+ a)) 

(sin 62 — sin 61 ) 



(cos 6*i sin 2 — sin 81 cos 9 2 ) 



(206) 



and 



(a) = 2C p (a) 



d9. 



y max J — 



d8 2 (cos 2p (6' 1 - a) ± cos 2p (6»i + a)) 
.ax — #1 ) ( cos 01 sin 02 — sin 6*i cos 9 2 ) 



sin 6*2 cos 2 0i — sin 0i cos 2 6*2 
sm.0! sin 2 (sin 0! — sin# 2 ) 
sin 6 1 sin 62 (cos 0i — cos 02 ) 



(207) 



The choice of signs indicated on the r.h.s. is to be understood as follows. The + sign is used for 
io, i\, and the — sign for i 2 . Using these expression all constants c\, C2, and Bq through B§ can 
be determined. 
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A. 5. Numerical values of the different coefficients 



71 
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< 1(T 9 


1.05321 


1.64421 


2.26475 


2.58598 


2.7831 


a 


0.185067 


0.250961 


0.374958 


0.594681 


0.751558 


0.862493 


6 


0.185067 


0.297586 


0.45714 


0.727492 


0.916382 


1.04853 


c 


0.432281 


0.308721 


0.315767 


0.368355 


0.416152 


0.452717 


c 


0.432281 


0.971317 


1.48443 


2.25966 


2.76618 


3.10848 


d 


-0.247214 


-0.246906 


-0.270197 


-0.311477 


-0.341938 


-0.364713 


r 


1.0 


0.914026 


0.438156 


-0.042155 


-0.260547 


-0.383077 




1.0 


0.843324 


0.820227 


0.81744 


0.820136 


0.822574 



TABLE I: The microscopic constants Co, ci and B M , fx = 0, . . . , 5, the entries a, b, c,c' ,d of the matrix Af 
calculated from the microscopic model for scattering for different bias intensities p where the maximum 
scattering angle is 6 max = 7r/2 — 0.01. 
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B 4 
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0.915888 


b 


0.339085 


0.501589 


0.681264 


0.952482 


1.1194 


1.23675 


c 


0.712094 


0.521519 


0.513996 


0.548739 


0.580671 


0.606106 


t 

c 


0.712094 


1.36669 


1.89275 


2.61839 


3.04756 


3.3414 


d 


-0.373009 


-0.370911 


-0.38917 


-0.419075 


-0.439207 


-0.453324 


r 


1.0 


0.868488 


0.57427 


0.252696 


0.0920036 


-0.00397738 


t 


1.0 


0.7989 


0.75906 


0.74107 


0.740293 


0.740561 



TABLE II: The microscopic constants co, ci and B^, y, = 0, . . . , 5, the entries a, b, c, c', d of the matrix A-f 
calculated from the microscopic model for scattering for different bias intensities p where the maximum 
scattering angle is 9 mal = n/2 — 0.05. 
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p 





1 


2 


4 


6 


8 


Cl - 1 


0.154395 


0.224047 


a ati /i oa 

0.271432 


0.32539 


0.355728 


0.37555 


C2 


< 10 


a ncooo /in 

-0.0528,349 


-0.0821605 


a 110cm 

-0.113501 


A 1 O A C C\C 

-0.130596 


-0.14166 


-Bo 


-0.199655 


a a/^tta a 

-0.267729 


a 1 1 tat 

-0.311707 


a a a r- 

0.360256 


-0.386527 


A A AO AO O 

-0.403233 


Bi 


1.79315 


1.70859 


1.69527 


1 f Ci TO V 

1.68785 


1 f Ci -1 

1.68312 


1 /""TATA 

1.67972 


T) 

B 2 


< 10 


A H000070 


a ao c a c 1 

0.0350513 


0.0396921 




A A A A 1 fl07 

U.U4^1(J27 


B 3 


< 10 


-0.0272506 


-0.0938402 


-0.16046 


A "I A A O AA 

-0.192802 


A A 1 1 O T A 

-0.211874 


T) 


< 10 


r\ A /I OAT /I 1 


a nooi ofi/i 

-0.0321304 


a fifinon 01 a 

0.000361819 


A AAA£JAAO 

0.0226003 


A AO O AO 1 

0.038231 


T) 

B5 


< 10 


0.0590451 


0.0753626 


A AOCO/1CO 

0.0858452 


A AOAOA^/1 

0.0893064 


A AAAOAOI 

0.0908931 


a 


0.224/0 


4.46981 


3.99399 


o.oooyo 


O 00(171 

3.334/ 1 


3.206/ / 


b 


5.22475 


5.48863 


5.38669 


5.23452 


5.14104 


5.08701 


c 


7.72907 


6.02669 


5.24234 


4.56791 


4.25716 


4.07678 


d 


7.72907 


8.43526 


8.55002 


8.60126 


8.61323 


8.63027 


d 


-2.50432 


-2.39596 


-2.11555 


-1.82208 


-1.68225 


-1.60082 


r 


1.0 


0.682741 


0.765062 


0.912647 


1.00577 


1.06836 


t 


1.0 


0.814376 


0.741456 


0.678371 


0.648644 


0.630384 



TABLE III: The microscopic constants Co, ci and fj. = 0, . . . , 5, the entries a, fo, c, c', d of the matrix A-j 
calculated from the microscopic model for scattering for different bias intensities p where the maximum 
scattering angle is 8 max = it/ 4. 



V 





1 


2 


4 


6 


8 


ci - 1 


0.0335067 


0.0428648 


0.0524673 


0.0641093 


0.0705982 


0.0747538 


C2 


< lO" 9 


-0.00668779 


-0.0121806 


-0.018075 


-0.020776 


-0.0221954 


Bo 


-0.0485058 


-0.0601859 


-0.0709311 


-0.0839805 


-0.0912928 


-0.0959171 


Bi 


1.94764 


1.89848 


1.86723 


1.86271 


1.86831 


1.87199 


B 2 


< lO" 9 


0.00497716 


0.00718661 


0.00839984 


0.00836975 


0.00806044 


B 3 


< lO" 9 


0.0189319 


-0.0293591 


-0.104785 


-0.149741 


-0.177526 


Bi 


<10" 9 


-0.0118706 


-0.0131302 


-0.00748694 


-0.000990488 


0.00438443 


B 5 


< lO" 9 


0.010374 


0.0158534 


0.0190278 


0.0189925 


0.018225 


a 


24.6907 


22.0583 


19.3127 


16.6004 


15.1812 


14.279 


b 


24.6907 


25.1971 


24.085 


22.3923 


21.0886 


20.0854 


c 


34.9988 


29.4735 


25.2402 


21.4431 


19.66 


18.5982 


c 1 


34.9988 


35.9891 


35.1548 


33.5005 


31.975 


30.7187 


d 


-10.308 


-10.0378 


-9.07808 


-7.69791 


-6.91559 


-6.43566 


r 


1.0 


0.697308 


0.677047 


0.784092 


0.890948 


0.973368 


t 


1.0 


0.87543 


0.801857 


0.741342 


0.719877 


0.710911 



TABLE IV: The microscopic constants Co, ci and £? M , (j, = 0, . . . , 5, the entries a, b, c,c' ,d of the matrix Af 
calculated from the microscopic model for scattering for different bias intensities p where the maximum 
scattering angle is 9 m ax = tt/S. 
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Appendix B: Response functions 
B.l. Region I 



The <?ij can be expressed as: 



r+oo r-too 

/ dq [a* 3 e- iqx + ai e tqx ]e lXiqz + dq [ale~ tqx + a 3 e iqx ]e iXsqz , (208) 
Jo Jo 

r+oo r+oo 

/ dq [(Xla 3 )*e- iqx + Xla A e %qx ]e %Xiqz + dq [(X 4 2 a 4 )* e~ lqx + X 2 a 3 e lqx ]e tX ' qz , 
Jo Jo 



r+oo r+oo 

- I dq[(X 3 a 3 )*er lqx + X 4 a 4 e iqx ]e iXiqz - / dq [(X 4 a A )*e- iqx + X 3 a 3 e lqx ]e iX3qz 



o 



o 



The top conditions H32I33[) allow to calculate the coefficients a 3 and 04. They read: 

f F 



"4 



X^ — X 3 2tt 



(X4 cos 6*o + sin^o), 



1 F 

— (X 3 cos0r, + sin0 o ) 



X 3 — X4 2tt 

To perform the integrals over q, it is useful to define the two following integrals: 

az =F i(3z 



+00 



I± = dq cos(gx) e- aqz±il3qz 



r+oo 

J± = dq sin(gx) e - aqz±l(3qz 
Jo 



(az =F i(3z) 2 + x 2 
x 

(az =F i(3z) 2 + x 2 



We then get 



2ir 2/3 
2tt 2(3 



/3 cos 8q + — h a cos 8q + - . h sin 9q - + 



2i 



2i 



(5(a 2 + 1 ) cos6> /+ +/ - a(a 2 + (3 2 )cos9 - + 



2i 



-(a 2 -[3 2 ) sin 9 a J+ n J - + 2a(3 sin 0„ J+ ± — 



2tt 2/3 



2i 

(a 2 + /3 2 ) cos O ^+-^ + /3 sin d 1 -±±^- 
Zl 2 



a sin Oq - 



I+-I- 



2i 



(209) 
(210) 

(211) 
(212) 

(213) 
(214) 

(215) 

(216) 
(217) 



B.2. Region II 



The <Jij can be expressed as: 
f +00 



r+oo r+oo 

/ dq \a\eT iqx + ai e lqx ]e lXiqz + dq[a 
Jo Jo 



<r M = / ^ [(X|a 4 )*e- l9a; + X 2 ai e lqx ]e tX * qz + / rfg [(X 2 a 3 )*e~ lqx + X 2 a 3 e* 9a V X39Z , 
Jo Jo 



(218) 



(219) 



r+oo r+oo 

/ dq [(X ia4 )*e~ lqx + X iai e lqx ]e tXiqz - dq [(X 3 a 3 )*e~ tqx + X 3 a 3 e lqx ]e tX3qz 
Jo Jo 



(220) 
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The top conditions i|32l33[) give again 



a 3 = „ 1 v ^(Xtcosflo + sinflo), (221) 

A 4 — A3 Z7T 

fl4 = v 1 y ^(A- 3 cosg + sing ). (222) 

A3 — A4 Z7T 



In this case, the useful integrals are 

r> + 00 



1(a) 



f + az 

/ dq co S (qx)e- a ^ = - 2 , (223) 

J (az) 2 + x l 



x 



We then get 



F 


2 


2tt ai 


— a\ 


F 


2 


2tt ai 


— a\ 


+a\a\ 


cos 


F 


2 



J(a) = / dq sm(qx)e- a " z = (224) 

(az) + x 



[a2 cos0o^(ai) + sin6>o J(«i) — ai cosflo-f (0:2) — sin6*o J(«2)] , 
\_-a\a2 cos 9 I (ai) — a\ sin9 J(ai)) 

il(a 2 ) + a 2 2 sin 6 J(a 2 )] , (225) 
[«i«2 cos #0 J(ct!i) — sin 6*o/(ai)) 



27T «2 — a i 

—a 2 a\ cos 9 J(a- 2 ) + 0:2 sin 6*0/(02))] . (226) 



35 



P.-G. de Gennes, Physica A 261, 293 (1998). 

S.B. Savage, in Physics of Dry Granular Media, H.J. Herrmann, J. P. Hovi and S. Luding, Eds., NATO 
ASI, 25 (1997). 

J.-P. Bouchaud, P. Claudin, M.E. Cates, J. P. Wittmer, in Physics of Dry Granular Media, H.J. 
Herrmann, J. P. Hovi and S. Luding, Eds., NATO ASI, 97 (1997). 

M.E. Cates, J. P. Wittmer, J.-P. Bouchaud and P. Claudin, Phil. Trans. Roy. Soc. Lond. A 356, 2535 
(1998). 

F. Cantelaube, J. Goddard, in Physics of Dry Granular Media, H.J. Herrmann, J. P. Hovi and S. 
Luding, Eds., NATO ASI, 123 (1997). 

P. Claudin, J.-P. Bouchaud, M.E. Cates and J. P. Wittmer, Phys. Rev. E 57, 4441 (1998). 

J.-P. Bouchaud, M.E. Cates, and P. Claudin, J. Phys. (France) I 5, 639 (1995). 

S. Ouaguenouni, J.N. Roux, Europhys. Lett. 32, 449 (1995); Europhys. Lett. 39, 117 (1997). 

C. Moukarzel, Phys. Rev. Lett. 81, 1634 (1998). 

V. Tkachenko and TA. Witten, Phys. Rev. E 60, 687 (1999) 

S.F. Edwards and D.V. Grinev, Phys. Rev. Lett 82, 5397 (1999); S.F. Edwards, D. Grinev, Physica 
A 294, 57 (2001). 

R.C. Ball and R. Blumenfeld, Phys. Rev. Lett. 88, 115505 (2002). 

J. P. Bouchaud, P. Claudin, D. Levine, and M. Otto, Eur. Phys. J E 4, 451 (2001). 

J.E.S. Socolar, P. Claudin, and D.G. Schaeffer, Eur. Phys. J. E 7, 353 (2002), and Erratum: Eur. 

Phys. J. E 8, 453 (2002). 

E. Clement, G. Reydellet, L. Vanel, D.W. Howell, J. Geng and R.P. Behringer, XBIth Int. Cong, on 
Rheology, Cambridge (UK), vol. 2, 426 (2000). 

G. Reydellet and E. Clement, Phys. Rev. Lett. 86, 3308 (2001). 

D. Serero, G. Reydellet, P. Claudin, E. Clement and D. Levine, Eur. Phys. J E 6, 169-179 (2001). 
C. Goldenberg and I. Goldhirsch, Phys. Rev. Lett. 89, 084302 (2002). 

C. Gay, R. da Silveira, cond-mat/0208155. R. da Silveira, G. Vidalenc and C. Gay, cond-mat/0208214. 
J. P. Wittmer, A. Tanguy, J.-L. Barrat, L. Lewis, Europhys. Lett. 57, 423 (2002). A. Tanguy, J. P. 
Wittmer, F. Leonforte, J.-L. Barrat, Phys. Rev. B 66, 174205 (2002). 
L. Landau and E. Lifshitz, Elasticity theory, Pergamon, New York (1986). 

C. Truesdell and W. Noll, in Handbuch der Physik III/3, S. Fliigge/C. Truesdell Eds., 1-579 (1965). 
C. Fichera, in Handbuch der Physik VIa/2, S. Fliigge/C. Truesdell Eds., 391 (1972). 
A.E. Green and G.I. Taylor, Proc. Roy. Soc. Lond. A 173, 162 (1939). A.E. Green, Proc. Roy. Soc. 
Lond. A 173, 173 (1939). 

L. Breton, P. Claudin, E. Clement and J.-D. Zucker, Europhys. Lett. 60, 813 (2002). 

J. P. Wittmer, P. Claudin, M.E. Cates and J.-P. Bouchaud, Nature 382, 336 (1996); J. P. Wittmer, P. 

Claudin, M.E. Cates, J. Phys. (France) I 7, 39 (1997). 

J. Geng, R. Reydellet, E. Clement, and R.P. Behringer, cond-mat/0211031, submitted to Physica D. 

E. Clement, R. Reydellet, B. Behringer, and J. Geng, private communication. Guillaume Reydellet, 
PhD Thesis, Mesure experimentale de la fonction de reponse d'un materiau granulaire, Universite 
Pierre et Marie Curie, Paris, (2002). 

J.E.S. Socolar, Discrete models of force chain networks, to appear in Discrete and Cont. Dyn. Systems 
Ser. B (2003). 

An alternative definition of P(f, n, r) such that P(f, n, r) ~ f~ D is used in [gsj which makes invariance 
of the full non-linear model with respect to rescaling of all force intensities most transparent. For the 
present work where mergings of force chains are neglected, this aspect is not relevant. 
M.C.W. van Rossum and Th.M. Nieuwenhuizen, Rev. Mod. Phys. 71, 313 (1999). 
J. P. Bouchaud, P. Claudin, M. Otto, J.E.S. Socolar, ongoing work. 

J. Gamier, PhD thesis, Tassement et contraintes. Influence de la rigidite de la fondation et de 
V anisotropic du massif., Universite de Grenoble (1973). 



